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The evolution of satellite galaxies is shaped by their constant interaction with the circum galactic medium 
surrounding central galaxies, which in turn may be affected by gas and energy ejected from the central 
supermassive black hole!~°. However, the nature of this coupling between black holes and galaxies is highly 
debated’~’ and observational evidence remains scarce!”!!, Here we report an analysis of archival data 
on 124,163 satellite galaxies in the potential wells of 29,631 dark matter halos with masses between 10! 
and 10'* solar masses. We find that quiescent satellites are relatively less frequent along the minor axis 
of their central galaxies. This observation might appear counterintuitive as black hole activity is expected 
to eject mass and energy preferentially in the direction of the minor axis of the host galaxy. However, we 
show that the observed signal results precisely from the ejective nature of black hole feedback in massive 
halos, as active galactic nuclei-powered outflows clear out the circumgalactic medium, reducing the ram 
pressure and thus preserving star formation. This interpretation is supported by the IllustrisTNG suite 
of cosmological numerical simulations, where a similar modulation is observed even though the sub-grid 
implementation of black hole feedback is effectively isotropic. Our results provide compelling observational 
evidence for the role of black holes in regulating galaxy evolution over spatial scales differing by several 
orders of magnitude. 


We use catalogs of groups and clusters!” to identify satellite galaxies in the Sloan Digital Sky Survey 
(SDSS)! data. For each satellite, its star formation rate (SFR) and stellar mass are measured!*:!° (see 
Methods for details). We proceed to characterize how satellites are distributed around their centrals in 
the plane of the sky. As shown in Fig. 1, the orientation of a given satellite galaxy can be defined as the 
difference between the photometric position angle of the central’s major axis and the position angle of 
the satellite with respect to the central. We retrieve the position angle of the central’s major axis from the 
SDSS imaging pipeline!® and we use sky coordinates to measure the position angle between centrals and 
satellites. Note that a galaxy does not need to be late type or of disky shape for a major axis to be defined, 
as in fact the great majority of the central massive galaxies in our SDSS sample are ellipticals. With this 
definition, a satellite located along the major axis of its central would have an orientation equal to 0° (or 
180°/360°) and, conversely, a satellite located along the minor axis would have an orientation of 90° (or 
270°). 

Given this metric for the angular distribution of satellites galaxies and the ancillary SFR and stellar 
mass measurements, Fig. 2 presents our main results. Panel (a) shows how the fraction of quiescent 
satellites depends on their orientation with respect to the central galaxy. Each point indicates the fraction 
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Figure 1. Orientation of satellite galaxies around centrals. For a given dark matter halo with a 
central galaxy, we define the orientation of each satellite as the difference between the position angle of 
the central’s major axis and the position angle of the satellite with respect to the central, as indicated in the 
figure. A satellite located along the major axis of the central would have an orientation of 0° (or 
180°/360°); conversely, a satellite located along the minor axis would have an orientation of 90° (or 270°). 


of quiescent galaxies among all satellites observed in a given orientation bin, regardless of the properties 
of the hosting halo. In practice this means that every group and cluster in our sample can contribute to 
each data point, depending on the location of its satellites. It is evident that the fraction of quiescent 
satellites is maximal along the major axis of the central galaxy and minimal along the minor axis. In 
order to quantitatively assess the significance of this modulation, we fit the observed signal with a cosine 
function with three free parameters: the median quiescent fraction, the amplitude of the modulation, 
and a re-scaling of the assumed error to account for any source of uncertainty. The modulation is well 
represented by a cosine function with an amplitude of 0.025+0.001 on top of a 0.421+0.001 average 
quiescent fraction. Complementary, panel (b) shows the iso-quiescent fraction contour (f,=0.42, the 
average of our sample) as a function of cluster-centric distance (normalized to the virial radius) and 
orientation angle. Along the major axis of the central galaxy, satellites are preferentially quenched at larger 
cluster-centric distances, while satellites located in the direction of the minor axis survive in relatively 
larger numbers as star-forming objects to much closer distances. Note that both panels in Fig. 2 are 
obtained by stacking every group and cluster in our sample into a single pseudo-cluster, where the major 
axis of every central galaxy is aligned in the same direction. Note also that, in principle, the signal should 
be fully symmetrical around the 0°-90°range. 

We have tested that this anisotropic modulation in the fraction of quiescent satellites, is statistically 
robust and exhibits a number of interesting features. As noted above the significance of a positive amplitude 
is well beyond 30 (0.025 + 0.001). Moreover, if we randomize the position angle of each central galaxy 
and we measure again the amplitude of the signal, we find in this case no modulation at all. At the same 
time, the signal is also robust against the assumed functional form for the surface brightness distribution 
while measuring the position angle of the central galaxy. The results shown in Fig. 2 are thus robust and 
truly dependent on the orientation of satellites. We also find that the amplitude of the signal increases for 
more massive centrals and for less massive satellites. Moreover, the signal becomes stronger for satellites 
closer to the center of their host halo. Interestingly, the signal also depends on the mass of the black hole 
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Figure 2. Anisotropic distribution of quiescent satellite galaxies in SDSS. Panel (a) shows how the 
fraction of quiescent galaxies depends on the orientation with respect to the central’s major axis. 
Quiescent galaxies are relatively enhanced along the direction of the major axis and conversely, relatively 
deficient in the direction of the minor axis. The signal is well represented by a cosine function with an 
amplitude of 0.025 (black solid line), and the shaded area represents the 1o confidence interval (40.001). 
Panel (b) represents the cluster-centric distance (normalized to the virial radius) at which the quiescent 
fraction equates 0.42 (the average across the adopted SDSS sample), as a function of orientation. The 
shape of the iso-quiescent contour demonstrates that satellite galaxies along the major axis of the central 
are quenched further out from the center of the cluster than those satellites located along the minor axis. 


in the center of a halo: at fixed halo mass, the amplitude of the observed signal is stronger if the black 
hole hosted by the central galaxy is more massive. All these additional tests are presented in the Methods 
section. 

To investigate the origin of the modulation, we make use of the IllustrisTNG suite of cosmological 
numerical simulations!’, where we find a similar behavior in the fraction of quiescent satellites. Fig. 3 is 
based on the outcome of the TNG100 run, using high-realism synthetic SDSS-like images!® to measure 
the centrals’ position angle. As in the SDSS data, the signal can be modeled by a cosine function with an 
amplitude of 0.032+0.004. Note that the selection function of SDSS satellites suffer from observational 
biases and completeness issues due to, e.g. fiber collisions. Hence, the absolute fraction of quiescent 
fractions might differ between SDSS and TNG100!°, particularly in groups and clusters. Therefore, in 
Fig. 3 we have subtracted the average quiescent fraction of both datasets. Furthermore, we have verified 
that the signal is in place in IllustrisTNG irrespective of whether the photometric position angle or the 
intrinsic stellar angular momentum of the centrals are used to identify their major axes. 

Two classes of distinct evolutionary mechanisms may be at the origin of the anisotropic distribution of 
quiescent satellites around central galaxies that we observe both in the SDSS data and in the IlustrisTNG 
cosmological numerical simulation. 

First, this could be the manifestation of a large-scale structure phenomenon, whereby satellites in 
groups and clusters form and evolve under different conditions than those in the field even before falling 
into the primary halos and thus, their properties and distribution can be affected by processes unrelated to 
their current group/cluster environment”’. 

Second, this could be the result of a (host) halo phenomenon, i.e. emerging because of the very 
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Figure 3. SDSS vs IllustrisTNG. Black symbols and lines are the SDSS observed data as in Fig. 2. 
Green symbols show the anisotropic distribution of quiescent satellites as measured in the IllustrisTNG 
100Mpc volume cosmological numerical simulation. In the latter the amplitude (0.032) is similar to that 
observed in SDSS, and the green shaded area indicates the 1o confidence interval (0.004). Since the 
absolute fraction of quiescent galaxies may differ between observed and simulated data, the mean 
quiescent fraction is subtracted to both datasets. Due to the lower number of satellites, [lustrisTNG 
measurements correspond only to the 0°-180°interval. To show them along with the SDSS data we have 
assumed that the 180°-360° behaves exactly as from 0°-180°one. 


interaction between satellite galaxies and their host halo. This could be for example due to energetic 
feedback from e.g. the super-massive black holes that reside at the center of massive halos, which in turn 
can have an impact beyond the central galaxy itself, modulating also the evolution of the surrounding 
satellites”! through the alteration of the physical properties of the intra cluster/group medium. 

In order to tests these two scenarios, we compare the signal measured in IllustrisTNG and shown in 
Fig. 3 with the signal measured from the previous generation, Illustris cosmological simulation?~. For the 
purposes of this work, the most noticeable difference between these two simulations is that IllustrisTNG 
includes an improved treatment of active galactic nuclei feedback, in particular in the low accretion rate 
regime**>. Thus, by comparing the two simulations, we are controlling for large-scale structure effects, 
isolating the role of black hole feedback in shaping the observed signal. We find that the modulation in 
the quiescent fraction significantly differs between the two simulations, being larger in IllustrisTNG (see 
Methods). For the first Illustris simulation, results are actually consistent with no modulation at a ~ 20 
level (0.0130.07). Furthermore, within the IllustrisTNG galaxy population, the anisotropic signal is 
driven by satellites that quenched within their current host halo (i.e. were not quenched prior to infall due 
to, e.g., pre-processing) and by the star-forming population of satellite galaxies (see Methods). 

We therefore conclude that the anisotropic behavior in the distribution of quiescent galaxies is a host 
halo phenomenon and, in particular, an observable manifestation of AGN feedback acting far beyond 
the extension of the central host galaxy**. This AGN-driven origin is particularly supported by the 
reported connection between the amplitude of the signal and the mass of the central black hole in the 
observed data, and more explicitly by the relation between energy injection and signal modulation in 
the simulations. Additionally, in both SDSS and IllustrisTNG data, the anisotropic quenching signal is 
stronger around quiescent rather than star-forming centrals, further suggesting a connection to the activity 
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of the super-massive black holes at the centers. 

Moreover, it is consistent with our analysis of the IllustrisTNG sample where, in the explored halo 
mass range, most satellite galaxies become quiescent after infalling in their current halo and the anisotropic 
signal in the simulation is dominated by such galaxies. Furthermore, AGN feedback is expected to impact 
more prominently the halo gas in the vicinity of the center, and in the SDSS data we find that satellites at 
smaller galactocentric distances exhibit a larger anisotropic modulation of quiescence. 
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Figure 4. Anisotropic CGM density in IllustrisTNG. Panel (a) illustrates how the mean relative gas 
overdensity changes as a function of orientation for central galaxies in the TNG100 simulation with stellar 
masses around the median of our SDSS sample (Meen ~ 10!! Mo). Averaged within the virial radius, the 
relative gas overdensity exhibits a similar behaviour as that measured for the fraction of quiescent galaxies 
in Fig. 2, as the circumgalactic medium gas density is relatively lower along the minor axis of centrals. 
Panel (b) shows a stacked average image of the gas overdensity (relative to its azimuthal average) over the 
same mass range, with the iso-quiescent fraction contour of Fig, | overlaid in white. For reference, the 
size of the virial radius is also indicated (dashed gray circle). 


The relative lower fraction of quiescent satellites along the minor axis of central galaxies and its 
relation with AGN feedback may at first seem counterintuitive because the energy and mass radiated 
by AGN activity are expected to escape the central galaxy preferentially along that direction. However, 
we argue that it is precisely the ejection of energy along the minor axis what drives the observed signal, 
as AGN feedback carves low-density bubbles in the CGM surrounding the central galaxy~>-*°. This is 
exemplified in Fig. 4, where we show how the average mass density of the CGM around central galaxies 
in IlustrisTNG displays exactly such anisotropic behavior, despite the fact that feedback processes are 
isotropic at the scale of energy injection. The interaction of the outflows with the galactic gas and the inner 
CGM is also likely responsible for the observed geometry, as the outflowing gas would tend to follow 
the path of least resistance, decoupled from the direction of energy injection at the smallest scales. This 
is supported by the fact that in the IlustrisTNG simulations, energy injection from super-massive black 
holes does not have a preferred direction, and yet outflows appear bi-polar on large scales. In any case, 
the anisotropic distribution of the CGM around central galaxies shown in Fig. 4, a direct IllustrisTNG 
prediction for the proposed scenario, should be observable in X-rays with sufficient statistics. 

All our findings support the idea that it is the interaction between satellites and the CGM, in turn 
modulated by AGN activity, that drives the observed signal. We suggest that, as satellite galaxies pass 
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through these low-density regions, processes directly responsible for their quenching such as ram pressure 
stripping become less efficient, increasing the relative abundance of star-forming galaxies along the 
direction of the minor axis. This interpretation is consistent with our analysis of the SDSS data, as ram 
pressure stripping is expected to affect more severely low mass satellites hosted by more massive centrals 
as indeed observed?’:*®. We cannot in fact exclude a different scenario, whereby the star formation activity 
of the satellites is enhanced rather than their quenching suppressed: it is possible that star-formation is 
enhanced along the outflowing material”, further increasing the fraction of star-forming satellites in the 
direction of the minor axis. 
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Methods 
1 Sample properties 


1.1 Sloan Digital Sky Server data 

Galaxy groups and clusters are selected!” from the 10th Data Release of the Sloan Digital Sky survey 1°. 
For each group and cluster, the Navarro-Frenk-White*” mass of the dark matter halo is provided, along 
with an estimation of the virial radius. In addition, the catalog also contains a list of satellites and identifies 
the central of each halo. We cross-match this catalog of groups and clusters with a catalog of SFR!° and 
stellar mass'* measurements based as well on spectro-photometric data from the SDSS. This matching of 
catalogs provides a final sample of 124,165 satellites with stellar masses in the M, ~ 10®?-10!!°M. range 
associated to 29,631 galaxy groups and clusters, hosted in halos with masses ranging from Mhalo ~ 1011-7 
to Mhalo ~ 10!45 Mo and with central galaxy mass in the M, ~ 10?-10!!°M. range, at a median redshift 
of z = 0.08. The typical distance of satellites with respect to central galaxies ranges from ~50 kpc to 
~ 900 kpc. No cuts were imposed on the basis of the b/a axis ratios of central galaxies. 

The separation between star-forming and quiescent satellites is based on their location with respect to 
the star formation main sequence (SFMS). In particular, we find that the SFMS in our sample of satellites 
is well-fitted*! by log SFR = 0.75logM, — 7.5. As described in the main text, given this definition of the 
SFMS and having stellar masses and SFR measurements for each satellite, we label as star-forming any 
satellite departing less than 1 dex from the SFMS. Conversely, a satellite is labeled as quiescent if its SFR 
is below the SFMS by more than 1 dex. 

In order to define the orientation of a satellite (see Fig. 1) the position angle (PA) of the central’s major 
axis has to be estimated. Thus, for each central in our sample, we use photometric information retrieved 
from the SDSS imaging pipeline!®. Since two measurements of the PA are given, depending on whether 
the central is better fit by an exponential or a de Vaucouleurs*” profile, we choose the parametrization 
with the highest likelihood given by the InLExp (exponential) or /nLDeV (de Vaucouleurs) keywords. To 
increase the robustness of the assumed quantities, we average likelihoods and PA over g, r, and i bands. 


1.2 IllustrisTNG simulation data 

We focus our analysis of the IlutrisTNG cosmological magneto-hydrodynamical simulations on the 100 
Mpc volume run, known as TNG100, for which synthetic SDSS-like images and photometric measurements 
have been made in a forward-modeling fashion!?. For simulated galaxies with stellar mass above 
M, = 10? Mo, we generate an SDSS-like synthetic image using the SKIRT radiative transfer code**** 
and perform a two-dimensional Sérsic fit using the statmorph code!8. We then use the best-fitting PA in 
the same way as in the SDSS data. This approach has two main advantages. First, image fitting is done 
on the projected X-Y plane of the simulation, which naturally emulates the distribution of orientations 
expected for central galaxies in the SDSS data. Second, a Sérsic fit generalizes the exponential and de 
Vaucouleurs profiles used in the SDSS imaging pipeline, ensuring a relatively fair comparison between 
observations and simulations. Images were synthesized using the snapshot 99 corresponding to redshift 
z = 0 and we adopt the same metric and definition of quiescent and star-forming as in the SDSS data: note 
that the SFMS of TNG100 is well consistent with observational constraints at low redshifts*>. 

Flagging out those centrals with unsuccessful Sérsic fits, we make use of a total of 880 halos (and 
centrals) in the host halo mass range of Mhalo ~ 10!* -10!4? and of 8,552 satellites with stellar masses 
ranging from M, ~ 108 to M, ~ 10116 Mo. At the resolution of TNG100, the less massive satellites 
would contain of the order of ~ 100 stellar particles!’. For central galaxies, stellar masses range from 
M, ~ 107° to 10!2 Mo. Details on galaxy and halo identification can be found in the IllustrisTNG 
presentation papers*!’. Despite the large cosmological volume of 100 Mpc, the number of objects is 
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significantly lower than in SDSS, although enough to reveal the presence of a modulation in the number 
of quiescent satellites (see Fig. 3). 

It is worth emphasizing that the goal of comparing SDSS and IllustrisTNG data is to understand the 
origin of the observed modulation in the fraction of quiescent galaxies, and not to provide an even-handed 
comparison between observed and simulated properties of galaxies in absolute terms: previous works 
have shown that the IllustrisTNG galaxy population is in good agreement with SDSS results, e.g. in terms 
of galaxy colors, global and small-scale stellar morphologies, SFRs, and quenched fractions! °°, In 
particular, the IlustrisTNG outcome is in striking quantitative agreement in comparison to SDSS data 
(differences smaller than < 5 — 10 percentage points) for global central and satellite quenched fractions at 
z= 0.1 and, in terms of satellites quenched fractions as a function of halocentric distance, for intermediate, 
group-mass scale hosts (10! — 10!4 M.)!°, which dominate the host mass distribution of both the SDSS 
and IllustrisTNG samples adopted here. We rely on those findings to elect IllustrisTNG as our simulation 
counterpart. At the same time, we note that an exact matching between the SDSS and IllustrisTNG 
samples is not critical for the purpose of this work, as we can get insights by focusing on the relative 
effects of the satellite locations and by marginalizing over the possibly-different absolute fractions of 
quiescent galaxies. On the other hand, matching exactly the galaxy samples would lead to an even lower 
number of available satellites, hampering the reliability of the comparison. 


2 Signal characterization 


In this section we detail the most meaningful tests done, both in the observed and simulated galaxy datasets, 
to understand the origin of the observed modulation in the number of quiescent satellites. 


2.1 Fitting the observed signal 

With no other motivation than quantifying the angular dependence of the fraction of quiescent galaxies, 
we fit the observed signal with a cosine function. In practice, we use a Bayesian MCMC sampler?’ to 
evaluate the following likelihood function 


In p(fq|0,4,b, f) = — +In(27s?) |, 
q l 


1 (fq, — a — bcos 26;)* 
ah |e 

where fg; is the observed fraction of quiescent galaxies at the 6; orientation, a is the median quiescent 
fraction, and b is the amplitude of the modulation. The error term is given by s? = o? + f*, where o is the 
estimated error, in our case estimated by bootstrapping, and f is the re-scaling term. 

This Bayesian framework allows us to explore the full posterior distribution and to naturally tests 
the null hypothesis (i.e. Zs the signal consistent with an amplitude equal to zero?). In Extended Data 
Fig. 1 we show the posterior distribution for the best-fitting solution of the SDSS data (Fig. 2). With these 
assumptions, we can reject the null hypothesis at a ~ 60 level. 

Note that error bars shown in all figures represent the combined s? = o? + f? uncertainty, and, for 
practical reasons, we fit for logarithmic re-scaling of the error In f. The inclusion of this additional error 
term f allows us to account for the various sources of error that might affect the observed data, such 
as the uncertainty on the estimated PA of the central galaxy, on the stellar mass and star formation rate 
measurements of the satellites, and on their stochastic distribution around centrals. 


2.2 Robustness of the signal against the assumed position angle 
Since the orientation of satellites strongly depend on the assumed PA of the central galaxy, we test 
the robustness of our results against errors and systematics on the PA determination. First, we test the 
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Extended Data Figure 1. SDSS posterior distributions for the best-fitting description of the 
angular modulation of satellite quiescence. We fit the observed data with a cosine function with three 
free parameters, the average quiescent fraction (f4), the amplitude of the modulation, and a re-scaling 
term for the expected error f. Posteriors are well-behaved and allowed us to reject the null-hypothesis at a 
~ 60 level. Blue solid vertical lines indicate the best-fitting values and the dashed ones the 10 confidence 


interval. 
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sensitivity of the observed signal to the functional form assumed for the surface brightness distribution 
of central galaxies. This is motivated by the fact that the light profile of galaxies is not accurately fitted 
by either a single exponential or de Vaucouleurs function as assumed by the SDSS photometric pipeline. 
Similarly to Fig. 2, panel (a) in Extended Data Fig. 2 shows the modulation in the fraction of quiescent 
satellites, but this time intentionally assuming the worst-fitting functional form according to SDSS. For 
example, if the SDSS photometric pipeline indicates that a galaxy is better fitted by an exponential profile, 
we selected the PA derived using a de Vaucouleurs, and vice versa. As it becomes obvious from panel (a) 
in Extended Data Fig. 2, the observed modulation in the fraction of quiescent satellites is insensitive to the 
functional formed used to fit the brightness profile of the central galaxy. 

Additionally, we also test if errors in the determination of the central’s PAs could have a significant 
impact on the observed signal. This is done by perturbing the PA of each central galaxy by a factor 6PA, 
drawn from a normal distribution ./ (0, APA”). By doing this we practically investigate the effect that a 
typical error of APA would have on the recovered signal, and it is shown in panel (b) of Extended Data 
Fig. 2. For errors on the individual PAs of up to 30°, there is still a clear modulation in the fraction of 
quiescent satellites. Note that the expected typical error in our sample of galaxies“? is of the order of ~ 2°, 
an order of magnitude smaller than the extreme test shown in Extended Data Fig. 2. Hence, we conclude 
that systematics and uncertainties related to the PA measurements based on the available photometric 
data in our sample of SDSS galaxies do not significantly affect our findings and thus our conclusions. 
For illustrative purposes, panels (c) and (d) in Extended Data Fig. 2 show examples of galaxies with de 
Vaucouleurs and an exponential light profiles, respectively. The uncertainty in the assumed position angle 
is indicated with the white dashed area. 

As a final test, we repeated our analysis but completely randomizing the PA of each central galaxy. 
In this case, the orientation of satellites becomes meaningless and therefore it is expected that the signal 
disappears. This is indeed the case as shown in panel (a) of Extended Data Fig. 3. Panel (b) in the same 
figure shows the posterior distribution, demonstrating that the amplitude of the signal vanishes when the 
centrals’ PAs are randomized. 


2.3 Dependence on galaxy properties in SDSS data 

The large number of objects in the SDSS sample of satellites allows us to investigate the dependence of the 
observed signal on different properties of centrals and satellite galaxies. Extended Data Fig. 4 summarizes 
the more meaningful trends we find in the SDSS data set. In particular, in panel (a) we split our sample 
into satellites close to the center of their halos/centrals (Rsat < 0.5Rxir) and satellites farther out from the 
center (Rgat > 0.5Ryir). For those satellites closer to the center we find an amplitude of 0.038 +0.002 in 
the modulation of the signal, while for those in the outskirts the modulation was weaker, with an amplitude 
of 0.023 +0.002. 

In addition, we also look at how the amplitude of the signal changes with the mass of the central 
galaxy, illustrated in Extended Data Fig. 4 panel (b). It is clear from this panel that the signal is stronger 
for satellites orbiting more massive centrals (log Mcen > 11 Mo), with an amplitude of 0.024 +0.002 
compared to what is observed for halos with a lower mass central (log Mcen < 11 Mo), where the amplitude 
is 0.007 +0.002. The opposite behavior is however observed when splitting our sample according to the 
mass of the satellites, as demonstrated in panel (c). Interestingly, the modulation of the signal is stronger 
for lower-mass satellites (log Msa < 10.5 Mo) with an amplitude of 0.027 +0.002. The amplitude 
measured for higher-mass satellites (log Msat > 10.5 Mo) is 0.014 +0.002. 

Here we also demonstrate what mentioned in the main text, namely that the strength of the observed 
modulation correlates with the (relative) mass of the super-massive black hole hosted by the central galaxy: 
see panel (d) of Extended Data Fig. 4. Following previous works*!, black hole masses are estimated 
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Extended Data Figure 2. Sensitivity of the SDSS signal to PA uncertainties. Panel (a) shows the 
fraction of SDSS quiescent galaxies as a function of the orientation based on the worst-fitting functional 
form (de Vaucouleurs vs. exponential) according to the SDSS photometric pipeline. The stability of the 
signal demonstrates that our results are robust against the photometric fitting procedure. In panel (b), 
colored curves indicate the best-fitting solution for SDSS data obtained while randomly perturbating the 
PA of the central galaxy by APA. For reference, black symbols and curves are the same as in Fig. 2. A 
clear modulation in the fraction of quiescent galaxies is observed even for APA~30°, which is an order of 
magnitude larger than the expected error on the individual PAs. Panels (c) and (d) show show the SDSS 
g-band images of galaxies best-fitted by a de Vaucouleurs (top row) and an exponential profile (bottom 
row), with the PA uncertainty indicated by the white shaded area. The adopted PA is indicated in the top 
left corner of each image. 


14/24 


Amplitude = 0.003*8-883 (b) 


@ Randomly-oriented centrals 


c 
[e] 0.001 
— fg = 0.4237; 
5 a 32388 
© o 
E s 
+ S 
j= ow 
© 7 
Y ko 
oO E 
5 4 D 
n o 
oO x logf = —4.046+9:957 
ci) 
5 >» 
iJ 
ž a 
> 
~ > 
xY : 
>” L— 
2S Sep P o © O a y Oo 49% 
0 50 100 150 _ 200 250 300 350 RR SEEE a 
Orientation [deg] A S amplitude o o po o logf 


Extended Data Figure 3. Test with randomized position angles. Panel (a) shows the fraction 
quiescent satellites in SDSS data after randomizing the PA of the central galaxies. As expected, no signal 
is recovered in this case. Panel (b) shows the posterior distributions for this test, where the modelled 
amplitude is consistent with no angular variation. 


using the empirical M,—o relation*’, and the stellar velocity dispersion of each central o is measured 
using the available SDSS spectroscopic data. At fixed halo mass, the amplitude of the observed signal is 
higher for those halos with a over-massive black hole in their centers (0.028 +0.002). For halos hosting 
under-massive black holes, the observed amplitude in the signal is 0.016 0.002. The expected uncertainty 
in these black hole mass estimates is of ~ 0.3 dex due the intrinsic scatter in the M.-o relation*!. In 
individual galaxies, this over-massive vs under-massive black hole metric has been now widely used to 
probe the interplay between black hole activity and star formation! !4!:-43-48, further supporting a black 
hole-related origin for the observed signal. For completeness, panels (e), (f), (g), and (h) in Extended Data 
Fig. 4 show the same analysis as in panels (a), (b), (c), and (d) but without removing the offset between 
the different sub-samples. 


2.4 Additional metrics for the satellite star formation properties 

Although we have mostly focused on the fraction of quiescent galaxies as a proxy for the properties of the 
satellite populations, alternatives measurements can be explored that could provide further insights on the 
origin of the signal. For example, the mean specific star formation rate at each orientation, or the average 
distance of satellites with respect to the star formation main sequence provide a less bimodal and more 
continuous characterization of the satellite population. The behavior of these two quantities is shown in 
Extended Data Fig. 5, and also exhibits the characteristic modulation reported for the fraction of quiescent 
galaxies. The consistency between this figure and the trend shown in Fig. 2 is an additional proof of the 
robustness of our results. 


2.5 Residual dependencies 

We have shown how the observed modulation in the fraction of quiescent satellites is stable and well- 
behaved. However, its peak-to-peak variation is only of ~ 5%, smaller than the expected change in the 
quiescent fraction when, for example, varying halo mass or cluster-centric distance. Thus, we also explored 
if the modulation could arise from variations in these properties. This was done by, first, characterizing the 
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Extended Data Figure 4. Characterization of the SDSS signal. In panel (a) we show that the 
modulation in the observed signal is higher for satellites closer to the centre (Rsat < 0.5Ryjir, orange 
symbols) than for those satellites in the outskirts (Rgat < 0.5Ryir, blue symbols). In panel (b) the signal is 
stronger for halos with more massive centrals (logMcen > 11Mọ, orange symbols) compared to the signal 
observed in halos with less massive centrals (logMcen < 11Mo, blue symbols). In panel (c) less massive 
satellites (logMsat < 10.5Mo, orange symbols) exhibit a larger variation than more massive ones 
(logMgat > 10.5Mo, blue symbols). Panel (d) reveals how the signal is also stronger in halos hosting more 
massive black holes in their center (orange symbols), compared to those with relatively less-massive 
central black holes (blue symbols) Panels panels (e), (f), (g), and (h) are equivalent to (a), (b), (c), and (d) 
but without removing the offset between the different sub-samples. 


16/24 


l 
> 
fon) 


a a 
x a 
n ki 
a 

pa 
2 cy 

oC 
= = 


log sSFR [y~+] 

I I | 
= o o 
o œ N 


| 
= 
H 


A Star formation main sequence [dex] 
b 
[Ce] 


—1.2 


0 50 100 150 200 250 300 £350 0 50 100 150 200 250 300 350 
Orientation [deg] Orientation [deg] 


Extended Data Figure 5. Alternative metrics for the characterization of SDSS satellites’ star 
formation status. The modulation observed in the average specific star formation rate (a) and distance 
from the star formation main sequence (b) of SDSS satellites closely follows that shown by the fraction of 
quiescent satellites in Fig. 2. Regardless of the metric used to characterize the star formation properties of 
satellite galaxies, there is a clear dependence on the orientation with respect to the central galaxy. Error 
bars indicate the 1o uncertainty and yellow lines mark the location of the minor and major axes. 


dependence of the quiescent fraction on Mhalo and Ryir with a quadratic polynomial fit. Then, we use these 
fits to estimate the amplitude of the signal that would result from a variation in halo mass or cluster-centric 
distance with orientation. The result of these tests is shown in Extended Data Fig. 6. 

For both halo mass and average radial distance we find subtle trends with orientation, but at a level 
much lower than it would be required to explain the modulation shown in Fig. 2. In particular, the average 
halo mass along the minor axis tends to be slightly higher (0.05 dex) than along the major axis. Since 
more massive halos tend to host more massive galaxies, the observed variation in halo mass which would 
lead to a variation in the fraction of quiescent satellites of 0.005+0.0001, much smaller than the observed 
one. Once corrected for halo mass, no other secondary trends are found. Regarding the typical radial 
distance, satellites along the minor axis are marginally farther away from the central that along the major 
axis. In this case the variation is even weaker than for the halo mass, with a best-fitting amplitude of only 
—0.0002 + 0.00005. 


2.6 Radial behavior 

Panel (b) in Fig. 2 demonstrates that the observed signal is radially well-behaved, as it allows to measure 
the contours of iso-quiescent fractions. For simplicity, Fig. 2 only shows one iso contour, but this idea 
can be further applied to different thresholds. In Extended Data Fig. 7 we show these contours at three 
different levels f, = {0.36,0.42,0.48} 


2.7 Illustris vs IllustrisTNG simulation comparison 

As an additional way to investigate the origin of the observed signal, Extended Data Fig. 8 illustrates 
how the outcomes of the IllustrisTNG and Illustris cosmological numerical simulations compare. To 
marginalize over the possible differences in the overall quenching of galaxies between Illustris TNG and 
Illustris*?, also here we subtract the average quiescent fraction of both datasets. 
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Extended Data Figure 6. Additional trends with halo mass and distance in SDSS. As in Fig. 2, 
black symbols represent the observed modulation on the SDSS data. The blue line indicates the change in 
the quiescent fraction that could be expected because of the average halo mass dependence on orientation, 
which is significantly smaller than the reported one. Similarly, satellites along the minor axis are 
marginally closer to the central galaxy than along the major axis, leading to a negative and even weaker 
modulation, as shown by the red line. 
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Extended Data Figure 7. Iso-quiescent fraction contours. Similarly to Fig. 4, the contours of 
constant fq are shown, but this time at three different levels (f4 = {0.36,0.42,0.48}). The background 
image corresponds to the IllustrisTNG gas over-density and the typical virial radius in the explored halo 
mass range is shown as a dashed grey circle. 
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As mentioned in the text, the TNG100 run of the IllustrisTNG series and the Illustris simulations both 
sample an approximately similar cosmological volume of 100 comoving Mpc: they in fact evolve the 
same initial conditions — but for small variations in the values of the adopted cosmological parameters — 
and have been performed at very similar numerical resolution. Illustris and IlustrisTNG differ in certain 
aspects of the underlying galaxy formation model’. In addition to an accurate treatment of the magneto- 
hydrodynamics within the simulation, IllustrisTNG improves upon the original Illustris galaxy formation 
model on three fronts (see Section 2.3 in?): chemical enrichment, galactic winds, and black hole feedback. 
The first one is unrelated to the topic of this work and the physical functioning of the stellar feedback 
under the form of galactic winds is similar albeit not identical in the two models. On the other hand, 
different mechanisms are invoked and implemented for the feedback from the super massive black holes 
in the two models, specifically at low accretion rates~. In particular, in both Illustris and IllustrisTNG, 
SMBH feedback is implemented by invoking three mechanisms? ?: thermal energy injection at high 
mass accretion rates, mechanical feedback at low mass accretion rates, and a sort of radiative feedback 
where gas cooling is further modulated by the radiation field of nearby AGNs. The two simulation 
models differ substantially only in the way the mechanical, low-accretion rate channel functions”: in 
IllustrisTNG, kinetic energy is injected into the surrounding gas as a pulsated wind, oriented in a different 
random direction at each SMBH timestep, and thus isotropic when averaged across any cosmological 
time scale of relevance; in contrast, in Illustris, thermal energy is injected in a highly bursty fashion into 
~ 50 — 100 kpc bubbles of gas at distances of a few 10 — 100s of kpc from the central galaxy. Such 
differences — whereby, in one case, AGN feedback affects gas at large radii rather than acting directly from 
the innermost regions of galaxies as in the case of IllustrisTNG -, result in important differences not only 
in the onset of galaxy quenching in central galaxies*>»*° but also in the properties of the circumgalactic 
medium (e.g. gas density)®?8-50-52, This implies that, in practice, by comparing in Extended Data Fig. 8 
Illustris and IlustrisTNG, we are effectively testing and varying the effects of different black hole feedback 
mechanisms on the satellite population. 

For the blue curve in Extended Data Fig. 8, Illustris galaxies are considered within similar ranges of 
host halo mass, and stellar mass of satellites and centrals as for the TNG100 sample (red curve). The 
signal is clear in IllustrisTNG, with an amplitude of 0.0320.004, whereas in Illustris the signal is weaker 
(0.013+0.007), consistent with almost no variation. Note that, since the number of satellites in these large 
simulated volumes is still small compared to SDSS, we symmetrize the data along 180°and the binning is 
also coarser. As the assembly of the larger scale-structure is the same in the two simulations, this piece of 
evidence favors the proposed scenario whereby feedback from super-massive black holes is responsible 
for the anisotropic quenching signal, through its effects on the properties of the circumgalactic medium. 

Note that since the mass function of Ilustris differs from that of Ilustris TNG®®, in principle, it could be 
possible that the larger amplitude in IllustrisTNG is due to a bias in the mass distribution even if galaxies 
were selected to cover the same mass range. Green symbols in Extended Data Fig. 8 show the variation in 
the fraction of quiescent IllustrisTNG satellites, but intentionally selected to reproduce the stellar mass 
distribution of satellites in Illustris. The fact that the amplitude of this signal remains unchanged for this 
M,-matched sample suggests that a mass bias is not responsible for the differences between Illustris and 
IllustrisTNG. 


2.8 Ejective AGN feedback at the origin of the quenching directionality in IllustrisTNG 

In IllustrisTNG, whether a massive galaxy is quenched or not is a direct indication of the effectiveness 
of black hole feedback at low accretion rates*°. Thus, a complementary way to asses the role of black 
hole feedback in shaping the observed modulation is by comparing the amplitude of the observed signal 
for quiescent and star forming centrals at fixed stellar mass. The outcome of this comparison is shown 
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Extended Data Figure 8. IllustrisTNG vs Illustris comparison. Modulation in the fraction of 
quiescent galaxies for the IlustrisTNG (namely, TNG100, red symbols) and the original Illustris (blue 
symbols) simulations. In green, the signal is shown for a sample of IllustrisTNG satellites with the same 
mass distribution as those in Illustris, to asses the possible effect of a mass bias between the two 
simulations. Both simulations probe a similar ~ 100 Mpc comoving cosmological volume and thus share 
the same large-scale structure properties, while the treatment of black hole growth and feedback is the 
most relevant difference between the two. However, it is clear that the amplitude of the modulation is 
significantly higher in IllutrisTNG (0.032+0.004) than in Illustris (0.0130.007). 


in Extended Data Fig. 9, panel (a), for central galaxies with stellar masses of log Mcen = 10.5 Mo. This 
particular mass range is selected so that both star-forming and quiescent centrals are abundant enough. 
Due to the additional constraints in stellar mass and SFR, the number of available satellites is rather small 
(typically ~ 100 per angular bin) and therefore the angular binning is coarser. 

It is apparent from Extended Data Fig. 9 that a systematic difference exists between the amplitude of 
the modulation in IllustrisTNG depending on the SFR of the central galaxy, as the signal is stronger for 
groups and clusters with a quiescent central galaxy, vanishing for satellites hosted by star-forming centrals. 
Since quiescence in massive IllustrisTNG galaxies is the consequence of an efficient and long-lasting 
AGN activity, the fact that the modulation is stronger for quiescent centrals further supports the black 
hole-related origin of the observed signal. Because of the low number of satellites in this test, however, 
the implications of Extended Data Fig. 9 should not be overstated. 

For completeness, we also investigate if the differences between quiescent and star-forming centrals 
shown in Extended Data Fig. 9, (a), is present in the SDSS observed data. As proved by Extended Data 
Fig. 9, panel (b), this in fact the case for SDSS central galaxies with stellar masses of log Mcen = 10.5 
Mo. 

In IllustrisTNG, the nature of the observed signal can be also tested by looking at how it depends on the 
amount of energy radiated, and on the mode in which this energy was radiated. In Extended Data Fig. 10, 
we select for TNG100 galaxies around log Meen = 10.5 Moand we split them in two groups, depending 
on the cumulative energy injected by their super-massive black holes. Panel (a) shows how the modulation 
in the quenching directionality signal is stronger for satellites around central galaxies whose black holes 
have injected, at a given stellar mass, more total (kinetic plus thermal) energy than the average: these 
results suggested by IllustrisTNG further link the observed signal with the black hole activity. Panel (b) 
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Extended Data Figure 9. Quiescent vs star-forming centrals in IllustrisTNG (a) and SDSS (b). 
In panel (a), at a fixed central stellar mass of ~ log Mcen = 10.5 Mo, the modulation in the fraction of 
quiescent satellites in TNG100 is shown for star-forming (blue) and quiescent (orange) centrals. Although 
with a limited number of satellites, the modulation in the signal appears to be stronger for quiescent 
centrals than for star-forming ones. Since quiescentness in IllustrisTNG is a strong indication of an 
effective black hole feedback, the fact that the signal is stronger for quiescent galaxies is also an indication 
of the proposed AGN-related origin for the observed quenching directionality. The modulation in the 
fraction of quiescent satellites is shown for star-forming (blue) and quiescent (orange) centrals in panel (b) 
but this time for SDSS galaxies, again of log Mcen = 10.5 Mo. The observed modulation is stronger for 
quiescent than for star-forming centrals as seen in IllustrisTNG. Solid lines and shaded areas indicate the 
best-fitting trends and 1-o confidence interval, respectively. 
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shows the result of splitting central galaxies but only according to the kinetic energy injected by their black 
holes when accreting at low rates. In this case, the signal almost vanishes for centrals that have undergone 
a relatively low amount of kinetic energy, once again supporting the idea that the signal is due to the 
ejective nature of black hole feedback, at least in the IlustrisTNG model. Whether, on the other hand, 
implementations of thermal and radiation SMBH feedback that are different from those adopted within 
Illustris and IllustrisTNG can reproduce the modulation of the fraction of quiescent satellites observed in 
SDSS remains to be determined. 

Complementary, panel (c) in Extended Data Fig. 10 shows the result of splitting IllustrisTNG galaxies 
by their relative black hole mass. As for the SDSS data (see panel (d) in Extended Data Fig. 4), the signal 
is stronger for satellites orbiting over-massive black holes centrals than for under-massive black hole 
ones, reinforcing the proposed connection between black hole activity and the observed modulation in the 
fraction of quiescent galaxies. 

In the observed SDSS data, the orientation angle between central and satellite galaxy is a projected 
one. That is, two satellites at apparently the same orientation may be actually located at a different 3D 
angle with respect to the central galaxy. Thus, the true underlying modulation is expected to be higher 
than the observed (projected) one. In order to explore this effect, we make use of the three-dimensional 
information provided by the IllutrisTNG simulation to de-project the observed modulation in the fraction 
of quiescent satellites. The result of this de-projection is shown in panel (d) of Extended Data Fig. 10. 
As expected from a three-dimensional effect such as the proposed black hole feedback mechanism, the 
amplitude of the signal, once de-projected, is larger than in the projected space. 

Finally, it is possible to use IllustrisTNG to control for the effect of large-scale structure in driving 
the observed quenching directionality by separating satellite galaxies according to when and where they 
quenched. In particular, as stated in the main text, the observed modulation in the fraction of quiescent 
satellites could in principle result from processes unrelated to their z ~ 0 host halo (e.g. pre-processing 
or quenching as centrals of a different halo). Satellite infall times have been studied for IllustrisTNG 
satellites’? and each satellite can be classified into four different groups: star-forming satellites, satellites 
quenched in their z ~ 0 host halo, pre-processed satellites (i.e. quenched as satellites in a different halo), 
and satellites quenched as centrals of a different halo (before being accreted into their z ~ 0 host). By 
construction, large-scale structure effects would be responsible for regulating the properties of the two last 
groups (i.e. pre-processed satellites and those quenched as centrals of a different halo), whereas black 
hole feedback from the central galaxy and the resulting quenching could only affect either satellites that 
quenched in their z ~ 0 host halo, as well as star-forming satellites (i.e. the first two groups). 

Panel (a) in Extended Data Fig. 11 shows the number of TNG100 satellites belonging to each group as 
a function of orientation. It is evident from this figure that, for the specific galaxy and host mass ranges 
adopted in this work*’, large-scale structure-sensitive satellites (red and orange symbols) are outnumbered 
by both star-forming, and particularly, satellites that quenched in their z ~ 0 host halo, and thus, large-scale 
effects could only affect a minority of the satellite population. 

Moreover, panel (b) in Extended Data Fig. 11 shows the modulation in the number of quiescent 
satellites but only taking into account star-forming satellites and those quenched within their z ~ 0 host 
halos, compared to the signal observed for the general population of IllustrisTNG satellites. After removing 
from the analysis those satellites that might be sensitive to large-scale structure effects, the strength of the 
signal remains the same, following the same trend represented in Fig. 3 and that observed in the SDSS data. 
We conclude therefore that the observed modulation in the fraction of IllustrisTNG quiescent satellites 
is indeed a host halo phenomenon, apparently related to the activity of supermassive black holes in the 
centers of groups and clusters. 
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Extended Data Figure 10. Dependence of the signal on the energy emitted by the super-massive 
black holes in IllustrisTNG. Panel (a) shows the fraction of quiescent satellites around centrals whose 
black holes have injected, relatively to their mass, more (red) and less (blue) total energy. Similarly, Panel 
(b) shows the same separation but in this case considering only the kinetic energy injected by the black 
holes. In both cases, the amplitude of the modulation is stronger when the total (panel a) and kinetic 
(panel b) energy released by the central black holes increase. Similar to Extended Data Fig. 4, panel (c) 
shows how the signal in IllustrisTNG depends on the relative mass of the central black hole, being 
stronger for more over-massive black hole galaxies. Finall, panel (d) shows, in red, the observed signal in 
IllustrisTNG, and in blue the de-projected using the underlying 3D satellite distribution. Note that in this 
panel we did not impose any cut in central stellar mass and therefore absolute values are different from the 
other panels. Error bars and shaded areas represent 1-o confidence intervals, and solid lines are the 
best-fitting solutions. 
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Extended Data Figure 11. Quenching directionality in IllustrisTNG. Panel (a) shows the number 
of TNG100 satellites in each orientation bin, depending on whether they are star-forming (blue symbols), 
quenched in their z ~ 0 host halo (green), were pre-processed and quenched in a different halo (orange), 
or quenched as centrals (red). The last two groups (red and orange symbols) are sensitive to large-scale 
structure effects, but only correspond to a small fraction of the total satellite population. In panel (b) the 
fraction of quiescent satellites as a function of orientation is shown but only for those satellites that 
quenched in their z ~ 0 host halo (green symbols). The amplitude of this modulation mimics that 
observed for all IllustrisTNG satellites (grey shaded area and black line). 
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